home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / ibeta.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  128 lines
  1. ;$Id: ibeta.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       IBETA
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the incomplete beta function, Ix(a, b).
  11. ;
  12. ; CATEGORY:
  13. ;       Special Functions.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;       Result = Ibeta(a, b, x)
  17. ;
  18. ; INPUTS:
  19. ;       A:    A positive scalar of type integer, float or double that 
  20. ;             specifies the parametric exponent of the integrand.
  21. ;
  22. ;       B:    A positive scalar of type integer, float or double that
  23. ;             specifies the parametric exponent of the integrand.
  24. ;
  25. ;       X:    A scalar, in the interval [0, 1], of type integer, float 
  26. ;             or double that specifies the upper limit of integration.
  27. ;
  28. ; EXAMPLE:
  29. ;       Compute the incomplete beta function for the corresponding elements
  30. ;       of A, B, and X.
  31. ;       Define the parametric exponents.
  32. ;         A = [0.5, 0.5, 1.0, 5.0, 10.0, 20.0]
  33. ;         B = [0.5, 0.5, 0.5, 5.0,  5.0, 10.0]
  34. ;       Define the the upper limits of integration.
  35. ;         X = [0.01, 0.1, 0.1, 0.5, 1.0, 0.8]
  36. ;       Allocate an array to store the results.
  37. ;         result = fltarr(n_elements(A))
  38. ;       Compute the incomplete beta functions.
  39. ;         for k = 0, n_elements(A)-1 do $
  40. ;           result(k) = Ibeta(A(k), B(k), X(k))
  41. ;       The result should be:
  42. ;         [0.0637686, 0.204833, 0.0513167, 0.500000, 1.00000, 0.950736]
  43. ;
  44. ; REFERENCE:
  45. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  46. ;       Cambridge University Press
  47. ;       ISBN 0-521-43108-5
  48. ;
  49. ; MODIFICATION HISTORY:
  50. ;       Written by:  GGS, RSI, September 1994
  51. ;                    IBETA is based on the routines: betacf.c, betai.c and
  52. ;                    gammln.c described in section 6.2 of Numerical Recipes,
  53. ;                    The Art of Scientific Computing (Second Edition), and is
  54. ;                    used by permission.
  55. ;-
  56.  
  57. function betacf, a, b, x
  58.   on_error, 2
  59.   eps   = 3.0e-7
  60.   fpmin = 1.0e-30
  61.   maxit = 100
  62.   qab = a + b
  63.   qap = a + 1.0
  64.   qam = a - 1.0
  65.     c = 1.0
  66.     d = 1.0 - qab * x / qap
  67.   if(abs(d) lt fpmin) then d = fpmin
  68.   d = 1.0 / d
  69.   h = d
  70.   for m = 1, maxit do begin
  71.     m2 = 2 * m
  72.     aa = m * (b - m) * x / ((qam + m2) * (a + m2))
  73.      d = 1.0 + aa*d
  74.      if(abs(d) lt fpmin) then d = fpmin
  75.      c = 1.0 + aa / c
  76.      if(abs(c) lt fpmin) then c = fpmin
  77.      d = 1.0 / d
  78.      h = h * d * c
  79.      aa = -(a + m) *(qab + m) * x/((a + m2) * (qap + m2))
  80.      d = 1.0 + aa * d
  81.      if(abs(d) lt fpmin) then d = fpmin
  82.      c = 1.0 + aa / c
  83.      if(abs(c) lt fpmin) then c = fpmin
  84.      d = 1.0 / d
  85.      del = d * c
  86.      h = h * del
  87.      if(abs(del - 1.0) lt eps) then return, h
  88.   endfor
  89.   message, 'Failed to converge within given parameters.'
  90. end
  91.  
  92. function gammln, xx
  93.   coff = [76.18009172947146d0,   -86.50532032941677d0,  $
  94.           24.01409824083091d0,    -1.231739572450155d0, $
  95.            0.1208650973866179d-2, -0.5395239384953d-5]
  96.   stp = 2.5066282746310005d0
  97.   x = xx
  98.   y = x
  99.   tmp = x + 5.5d0
  100.   tmp = (x + 0.5d0) * alog(tmp) - tmp
  101.   ser = 1.000000000190015d0
  102.   for j = 0, n_elements(coff)-1 do begin
  103.     y = y + 1.d0
  104.     ser = ser + coff[j] / y
  105.   endfor
  106.   return, tmp + alog(stp * ser / x)
  107. end
  108.  
  109. function ibeta, a, b, x
  110.   on_error, 2
  111.   if (x lt 0 or x gt 1) then message, $
  112.     'x must be in the interval [0, 1].'
  113.  
  114.   if (a le 0 or b le 0) then message, $
  115.     'a and b must be positive scalars.'
  116.  
  117.   if (x eq 0  or x eq 1) then bt = 0.0 $
  118.   else $
  119.     bt = exp(lngamma(a + b) - lngamma(a) - lngamma(b) + $
  120.              a * alog(x) + b * alog(1.0 - x))
  121.     ;bt = exp(gammln(a + b) - gammln(a) - gammln(b) + $
  122.     ;         a * alog(x) + b * alog(1.0 - x))
  123.   if(x lt (a + 1.0)/(a + b + 2.0)) then return, $
  124.     bt * betacf(a, b, x) / a $
  125.   else return, 1.0 - bt * betacf(b, a, 1.0-x) / b
  126. end
  127.  
  128.